/* ****************************************************************************
 *
 * Copyright (c) Microsoft Corporation. 
 *
 * This source code is subject to terms and conditions of the Microsoft Public License. A 
 * copy of the license can be found in the License.html file at the root of this distribution. If 
 * you cannot locate the  Microsoft Public License, please send an email to 
 * dlr@microsoft.com. By using this source code in any fashion, you are agreeing to be bound 
 * by the terms of the Microsoft Public License.
 *
 * You must not remove this notice, or any other, from this software.
 *
 * THIS CODE AND INFORMATION ARE PROVIDED "AS IS" WITHOUT WARRANTY OF ANY
 * KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE
 * IMPLIED WARRANTIES OF MERCHANTABILITY AND/OR FITNESS FOR A
 * PARTICULAR PURPOSE.
 * ***************************************************************************/


#if !defined(AFX_FIELDREGLATMAP_H__DDC727AC_69B9_4AEF_8EBF_422EDD9D0802__INCLUDED_)
#define AFX_FIELDREGLATMAP_H__DDC727AC_69B9_4AEF_8EBF_422EDD9D0802__INCLUDED_

#if _MSC_VER > 1000
#pragma once
#endif // _MSC_VER > 1000

#include "FieldMap.h"

namespace contour
{

/*! \brief A field objet on a regular latice

\ingroup FieldGroup

*/
template <class T>
class TFieldRegLatMap : public TFieldMap<T>  
{
public:
	//! \name Constructors
	//@{
	TFieldRegLatMap(const std::TSTLMatrix<T>& vData= std::TSTLMatrix<T>(0,0)):TFieldMap<T>(vData){};
	virtual ~TFieldRegLatMap(){};
	//@}

	//! \name Coordinates
	//@{
	//! Setting x coordinates
	void SetXCoord(const std::vector<T>& _vX);
	//! Setting y coordinates
	void SetYCoord(const std::vector<T>& _vY);
	//@}

	//! \name Eval function
	//@{
	//! Interpolate the map at x,y
	T Eval(T x, T y);
	//@}
private:
	int Find(int axe, T val) const;
	std::vector<T> m_vX;
	std::vector<T> m_vY;
	int BinSearch(const std::vector<T>& v, T val) const;
	double m_dVal;
};

template<class T>
void TFieldRegLatMap<T>::SetXCoord(const std::vector<T>& _vX)
{
	ASSERT(_vX.size()>= GetNCols());
	m_vX=_vX;
};

template<class T>
void TFieldRegLatMap<T>::SetYCoord(const std::vector<T>& _vY)
{
	ASSERT(_vY.size()>=GetNRows());
	m_vY=_vY;
};

template<class T>
T TFieldRegLatMap<T>::Eval(T x, T y)
{	
	int xi, yi;
	T alphaxi;
	T alphayi;

	// truncating x,y
	x=__max(m_vX[0], __min( m_vX[GetNCols()-1] , x) );
	y=__max(m_vY[0], __min( m_vY[GetNRows()-1], y));

	// getting index
	xi=Find(0,x);
	yi=Find(1,y);

	alphaxi=(x-m_vX[xi])/(m_vX[xi+1]-m_vX[xi]);
	alphayi=(y-m_vY[yi])/(m_vY[yi+1]-m_vY[yi]);

	if ( (xi == GetNCols()-1) && (yi == GetNRows()-1) )
	{
		return m_mMap(xi,yi);
	}
	else if ( xi == GetNCols()-1 )
	{
		// interpolating
		return m_mMap(xi,yi)*(1-alphayi)+m_mMap(xi,yi+1)*alphayi;
	}
	else if (yi == GetNRows()-1)
	{
		// interpolating
		return m_mMap(xi,yi)*(1-alphaxi)+m_mMap(xi+1,yi)*alphaxi;
	}
	else
	{
		// interpolating
		return m_mMap(xi,yi)*(1-alphaxi)*(1-alphayi)
			+m_mMap(xi,yi+1)*(1-alphaxi)*alphayi
			+m_mMap(xi+1,yi+1)*alphaxi*alphayi
			+m_mMap(xi+1,yi)*alphaxi*(1-alphayi);
	}
}

template<class T>
int TFieldRegLatMap<T>::Find(int axe, T val) const
{
	switch(axe)
	{
	case 0:
		return BinSearch(m_vX, val);
		break;
	case 1:
		return BinSearch(m_vY, val);
		break;
	default:
		ASSERT(false);
	}
	return 0;
};

template<class T>
int TFieldRegLatMap<T>::BinSearch(const std::vector<T> &v, T val) const
{
	int begin, end,mid;

	begin=0;
	end=v.size()-1;
	while ( !( (begin==end) || (begin==end-1) ) )
	{
		mid=begin+(end-begin)/2;
		
		if (val< v[mid])
		{
			end=mid;
		}
		else
		{
			begin=mid;
		}
	}
	return begin;
};

};

#endif // !defined(AFX_FIELDREGLATMAP_H__DDC727AC_69B9_4AEF_8EBF_422EDD9D0802__INCLUDED_)
